clc;clear;

% Load the MAT file
data = load('sph01_321_ref_par.mat');

% Access the matrix stored in Param_stack
params = data.Param_stack;

% Extract the 17th row (least square values)
least_squares = params(17, :);  % 17th row (since MATLAB is 1-based indexing)

% Sort the columns based on the 17th row values (least square values)
[~, sorted_indices] = sort(least_squares);

% Get the first 10 sorted columns
sorted_params = params(:, sorted_indices(1:20));

% Assign the first 10 sorted parameter sets to variables p1 to p10
p1 = sorted_params(:, 1);
p2 = sorted_params(:, 2);
p3 = sorted_params(:, 3);
p4 = sorted_params(:, 4);
p5 = sorted_params(:, 5);
p6 = sorted_params(:, 6);
p7 = sorted_params(:, 7);
p8 = sorted_params(:, 8);
p9 = sorted_params(:, 9);
p10 = sorted_params(:, 10);
p11 = sorted_params(:, 11);
p12 = sorted_params(:, 12);
p13 = sorted_params(:, 13);
p14 = sorted_params(:, 14);
p15 = sorted_params(:, 15);
p16 = sorted_params(:, 16);
p17 = sorted_params(:, 17);
p18 = sorted_params(:, 18);
p19 = sorted_params(:, 19);
p20 = sorted_params(:, 20);

% Initialize an array to store the A values
% Initialize an array to store the A values
A_values = zeros(20, 1);

% Define z range and C (assuming C is already defined)
z = (-200:1:50)';  % z range
C(2) = 1;  % Example, replace with your actual C value

% Perform the operation for p1
P = p1;  % Use p1 for the first calculation

% Modify the parameters for p1
Pnew = P;
Pnew(7) = P(4);  % Set Pnew(7) to P(4)
Pnew(10) = P(4); % Set Pnew(10) to P(4)
Pnew(13) = P(4); % Set Pnew(13) to P(4)

% Calculate ED profiles for p1
rho = fitref_ED_3box(P, C, z);       % Original ED profile for p1
rhonew = fitref_ED_3box(Pnew, C, z); % Modified ED profile for p1

% Calculate the difference and sum to get A for p1
A_p1 = sum(rho - rhonew);

% Display the sum for p1
disp(['Sum of differences (A) for p1: ', num2str(A_p1)]);

% Now loop through each parameter set (p1 to p10) and calculate A values
for i = 1:20
    % Get the current parameter set (p1 to p10)
    P = eval(['p', num2str(i)]);
    
    % Modify the parameters as per the instructions
    Pnew = P;
    Pnew(7) = P(4);  % Set Pnew(7) to P(4)
    Pnew(10) = P(4); % Set Pnew(10) to P(4)
    Pnew(13) = P(4); % Set Pnew(13) to P(4)
    
    % Calculate ED profiles using fitref_ED_3box
    rho = fitref_ED_3box(P, C, z);      % Original ED profile
    rhonew = fitref_ED_3box(Pnew, C, z); % Modified ED profile
    
    % Calculate the difference and sum to get A
    A = sum(rho - rhonew);
    disp(A);
    
    % Store the A value
    A_values(i) = A;
end

% Calculate the average and standard deviation of A values
A_avg = mean(A_values);
A_std = std(A_values);

% Display the results
disp(['Average of A values: ', num2str(A_avg)]);
disp(['Standard deviation of A values: ', num2str(A_std)]);

